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FOREWORD 


The  introduction  of  intership  data  networking  into  the  process  of  target  tracking  by 
multiple  platform  sensors  has  raised  the  question  of  how  to  deal  with  remote  sensor  data 
that  has  arrived  for  fusion  processing  after  a  delay  due  to  network  loading.  Conditions 
frequently  arise  wherein  standard  Kalman  filtering  formulations  must  be  modified  to  be  able 
to  accept  and  meaningfully  process  late  data.  This  report  presents  results  of  research  into 
the  relative  efficacy  of  a  complex,  rationally  derived  modified  Kalman  filter  update  versus 
faster  and  simpler  but  ad  hoc  Kalman  filter  modifications  for  accepting  late  data.  This  work 
was  conducted  under  the  sponsorship  of  the  Naval  Surface  Warfare  Center  Dahlgren  Division 
(NSWCDD)  AEGIS  Program  Office. 

This  document  has  been  reviewed  by  Dr.  T.R.  Rice,  Senior  Physicist,  Combat  Systems 
Technology  Group. 


Approved  by: 

DAVID  B.  COLBY,  Head 
Systems  Research  and  Technology 
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ABSTRACT 


When  techniques  for  target  tracking  are  expanded  to  make  use  of  multiple  sensors  in 
a  multiplatform  system,  the  possibility  of  time  delayed  data  becomes  a  reality.  When  a 
discrete-time  Kalman  filter  is  applied  and  some  of  the  data  entering  the  filter  are  delayed, 
proper  processing  of  these  late  data  is  a  necessity  for  obtaining  an  “optimal”  estimate  of  a 
target’s  state.  If  this  problem  is  not  given  special  care,  the  quality  of  the  state  estimates 
can  be  degraded  relative  to  that  quality  provided  by  a  single  sensor.  A  negative-time  update 
technique  is  developed  using  the  criteria  of  minimum  mean-square  error  (MMSE)  under  the 
constraint  that  only  the  results  of  the  most  recent  update  are  saved.  The  performance  of  the 
MMSE  technique  is  compared  to  that  of  the  ad  hoc  approach  employed  in  the  Cooperative 
Engagement  Capabilities  (CEC)  system  for  processing  data  from  multiple  platforms.  It 
was  discovered  that  the  MMSE  technique  is  a  stable  solution  to  the'  negative-time  update 
problem,  while  the  CEC  technique  was  found  to  be  less  than  desirable  when  used  with  filters 
designed  for  tracking  highly  maneuvering  targets  at  relatively  low  data  rates.  The  MMSE 
negative-time  update  technique  was  found  to  be  a  superior  alternative  to  the  existing  CEC 
negative-time  update  technique. 
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CHAPTER  1 
INTRODUCTION 


When  techniques  for  target  tracking  are  expanded  to  make  use  of  multiple  sensors  in 
a  multiplatform  system,  the  possibility  of  time  delayed  data  becomes  a  reality.  When  a 
discrete-time  Kalman  filter  is  applied  and  some  of  the  data  entering  the  filter  are  delayed, 
proper  processing  of  these  late  data  is  a  necessity  for  obtaining  an  “optimal”  estimate  of  a 
target’s  state.  If  this  problem,  which  may  arise  whenever  data  from  dispersed  sensors  are 
fused,  is  not  given  special  care,  the  quality  of  the  state  estimates  can  be  degraded  relative 
to  that  quality  provided  by  a  single  sensor. 

The  optimal  solution  to  this  problem  of  processing  late  data  requires  the  filter  algorithm 
to  return  to  the  state  estimate  and  error  covariance  that  occurred  just  prior  to  the  time 
corresponding  to  the  late  measurement.  Then  the  filter  process  is  restarted  from  the  time  of 
the  late  measurement.  Thus,  all  the  measurement  information  in  the  time  interval  between 
the  most  recent  measurement  and  the  late  measurement  must  be  saved  to  complete  this 
filtering  of  past  data.  In  order  to  bound  the  amount  of  information  to  be  saved,  a  maximum 
for  the  delay  time  must  be  declared.  Thus,  all  measurements,  state  estimates,  and  error 
covariances  within  the  delay  interval  must  be  saved.  This  filtering  approach  that  minimizes 
the  mean-square  error  may  not  be  practical  in  many  cases.  For  example,  when  a  tracking 
system  is  required  to  maintain  current  state  estimates  on  many  targets,  storing  the  past 
information  can  become  cumbersome. 

In  this  report,  the  problem  of  updating  the  state  estimate  with  late  data  is  considered 
for  the  case  where  only  the  most  recent  state  estimate  and  error  covariance  is  saved.  A 
solution  is  derived  for  this  problem  using  a  minimum  mean-square  error  (MMSE)  criteria. 
While  the  target  state  estimate  provided  by  this  approach  will  be  inferior  to  that  provided 
by  the  optimal  MMSE  approach  that  is  described  above,  this  suboptimal  approach  involves 
only  modifying  the  state  estimate  produced  with  the  most  recent  measurement  using  the 
late  measurement. 

Figure  1-1  illustrates  a  case  wherein  one  of  two  sensors  sends  data  to  the  site  of  fusion 
by  way  of  a  data  network.  In  the  illustrated  network,  the  input  queue  is  susceptible  to 
being  crowded  by  other  traffic,  causing  measurements  from  sensor  2  to  arrive  at  the  site 
of  data  fusion  late.  Due  to  the  late  arrival,  a  filter  update  is  made  using  measurements 
from  only  sensor  1,  and  if  the  measurement  from  sensor  2  is  to  be  utilized,  it  will  require 
special  treatment.  In  this  report,  the  process  of  using  late  measurements  to  improve  the 
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Figure  1-1.  A  Kalman  Filter  Receiving  Delayed  Measurements 

current  state  estimate  will  be  referred  to  as  a  negative-time  update;  and  the  time  interval 
between  the  most  recent  measurement  and  the  late  measurement  will  be  referred  to  as  the 
negative-time  interval.  Some  authors  refer  to  the  subject-matter  as  “data  senescence.” 

The  need  for  management  of  negative-time  intervals  arises  in  the  tracking  techniques 
used  for  implementation  of  Cooperative  Engagement  Capabilities  (CEC).  It  was  discovered 
in  the  Software  Requirements  Specification  (SRS)  document  of  CEC  [1],  that  an  ad  hoc 
approach  for  handling  negative-time  measurements  is  used.  A  major  motivation  for  under¬ 
taking  this  study  was  to  examine  the  efficiency  and  robustness  of  the  ad  hoc  treatment  of 
delayed  data  within  its  environment  -  an  extended  Kalman  filter  of  the  sort  described  in 
Appendix  A.  An  MMSE  approach  for  handling  negative-time  updates  is  derived  and  used  as 
a  baseline  for  evaluating  the  specified  CEC  approach.  The  specified  CEC  updating  method 
will  be  simulated  as  a  separate  program  and  will  not  be  tested  within  a  full  simulation  of 
the  CEC  system. 

This  report  is  organized  in  the  following  manner:  Chapter  2  discusses  the  problem  for¬ 
mulation  and  background  information.  Chapter  3  presents  the  development  and  derivation  of 
the  MMSE  negative-time  approach  along  with  a  simple  numerical  example  of  its  implemen¬ 
tation  and  an  analysis  of  constant  velocity  tracking  with  a  negative-time  update.  Chapter  4 
discloses  an  interpretation  of  the  CEC  approach  with  a  numerical  example  and  analysis  to 
parallel  that  in  Chapter  3.  Chapter  5  contains  details  and  results  of  Monte  Carlo  simulations 
used  to  evaluate  the  performance  of  the  ad  hoc  approach  in  a  simplified  but  realistic  context. 
Chapter  6  presents  conclusions  of  the  analysis. 
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CHAPTER  2 

BACKGROUND  AND  PROBLEM  FORMULATION 


In  this  study,  an  extended  Kalman  filter  is  used  to  update  a  multiple  source  track  file  by 
processing  measurements  of  range,  bearing,  and  elevation.  The  study  employs  a  simplified 
two-ship  communications  scenario  in  which  an  cwnship  sensor  is  tracking  a  maneuvering 
target  while  also  receiving  track  contacts  from  a  remote  ship  sensor  tracking  the  same  target. 
The  two  ships  are  confined  to  the  (East- West)  x-axis  in  ordtr  to  eliminate  the  need  for 
reference  frame  rotation;  and  the  earth’s  shape  is  considered  locally  flat  for  simplicity. 


Background 


In  CEC,  every  Cooperative  Unit  (CU)  platform  maintains  target  tracks  in  its  own  local 
coordinate  system.  Each  CU  uses  its  own  sensors  and  broadcasts  all  of  its  sensor  measure¬ 
ments  in  its  own  coordinate  system.  Periodically,  each  CU  broadcasts  information  on  its  own 
position.  When  ownship  receives  a  remote  sensor  measurement,  it  must  take  into  account 
the  coordinate  system  conversion  when  formulating  the  Kalman  filter  innovation  term  for 
tracking  at  ownship.  The  dynamics  and  measurement  models  assumed  for  the  target  in  track 
are  given  by 


Xk+i  =  F*A*  +  Gk  to*  (2.1) 

Zk  =  hk(Xk)  +  Vk  (2.2) 


where  to*  ~  7V(0,  Qk)  is  the  process  error  and  u*  ~  N(0,  Rk)  is  the  measurement  error,  re¬ 
spectively.  Matrix  F*  is  the  state  transition  matrix,  <7*  is  the  process  error  input  matrix,  and 
the  target  state  A*  contains  the  position,  velocity,  and  acceleration  of  the  target.  A  linear 
measurement  equation  can  be  defined  by  linearizing  the  nonlinear  observation  function.  The 
measurement  equation,  as  a  function  of  the  actual  target  state,  is  given  by 
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(2.3) 


x'k  =  Xfc-Xo* 

(2.4) 

y'k  -  Vk  -  Vob 

(2.5) 

Zk  =  2k  ~  Zob 

(2.6) 

and  (x/t,  z/c)  correspond  to  the  (x,y,  z)  coordinates  of  the  target  in  the  ownship  reference 
frame  at  time  k,  and  (x0{,,  y0(,,  zof))  correspond  to  the  origin  of  the  remote  ship  in  terms  of  the 
ownship  reference  frame  at  time  k.  Since  we  are  neglecting  motion  between  the  two  ships, 
x0b,  Vob,  and  z0 &  are  constants.  For  ownship  measurements,  (x06,  y0b,  z0b)  =  (0,0,0).  The 
r*,  bk,  and  e*  are  the  measured  range,  bearing,  and  elevation  relative  to  the  remote  sensor 
reference  frame  at  time  k.  The  vTk,  vk,  and  vk  are  the  measurement  errors  for  the  range, 
bearing,  and  elevation  at  time  k,  respectively. 

For  a  filter  processing  range,  bearing,  and  elevation  measurements, 


Using  a  first-order  Taylor  series  to  linearize  (2.7)  about  the  predicted  state  estimate  re¬ 
sults  in  the  linearized  measurement  matrix  given  by 
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where 


Hk  = 


_/ 

o 

Rk\k-i 

Rk\k-i 

0 

xk\k-l 

Rhk  |*-i 

Rhk\k-i 

zk\k-\x'k\k-\ 

0 

0  o 

0  0  0 


Xk  =  [  xk  xk  yk  yk  zk  zk  f 


Rk\k-i  =  y/K\k-i)2  +  (Vfclfc-i)2  +  (**|*-i)2 
Rhk\k-i  =  >/(**l*-»)a  +  (y*l*-i)a 


The  equations  for  the  extended  Kalman  filter  are  given  as  follows. 


(2.8) 


(2.9) 

(2.10) 

(2.11) 


Time  Update: 


=  P fc— i|fc— l  (2-12) 

Pk\k-1  =  Fk-lPk~l\k-lPk-l  +  Gk-lQk-lGk-i  (2.13) 

Measurement  Update: 

Xk\k  =  Xk\k-i  +  Kk[Zk  -  hk{Xk^)\  (2.14) 

Pk]k  =  [I  -  KkHk)Pk\k-i  (2.15) 

Kk  =  Py^HliHkP^HU  Rk\-'  (2.16) 


The  subscript  notation  (fc|j)  refers  to  the  estimate  at  time  k  given  measurements  through 
time  j.  The  Pk\k  is  the  error  covariance  associated  with  the  state  estimate  Xk\k,  and  Kk  is 
the  Kalman  gain. 


Problem  Formulation 


In  Figure  2-1,  a  time  line  is  given  for  a  situation  in  which  a  negative-time  update  of 
a  discrete  time  extended  Kalman  filter  is  necessary.  At  tk-2,  a  normal  filter  update  using 
timely  data  Zk...i,  results  in  state  estimate  Xk-2\k-2-  At  tk-i,  a  remote  sensor  obtains  a 
measurement  appends  these  data  with  the  relevant  time  tag  tk~i,  and  submits  a  data 
message  to  a  data  network  where  its  transmission  is  delayed.  At  tk,  another  measurement 
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time 


tk-  2 


tk- 1 


tk 


Figure  2-1.  Time  Line  for  Events  Requiring  a  Negative-Time  Update 

Zk  from  the  ownship  sensor  arrives  on  time  and  is  used  to  update  the  state  estimate  from 
Xk-2\k-2  to  Xk\k-  After  this  update  takes  place,  at  tk,  measurement  Zk-\  with  time  tag  f*_i 
arrives  after  its  delayed  transmission.  An  approach  to  incorporate  this  delayed  data  into 
an  “improved”  state  estimate  Xk\k+  is  developed  in  the  next  chapter.  A  related  problem 
has  been  studied  by  Thomopoulos  and  Zhang  [2],  wherein  data  reaches  a  fusion  site  after  a 
delay  whose  duration  is  known  statistically  but  not  exactly  (that  is,  when  no  time  tag 
is  appended  to  the  measurement).  Their  formulations  bear  some  resemblance  to  the  results 
obtained  in  the  next  chapter. 
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CHAPTER  3 

MINIMUM  MEAN-SQUARE  ERROR  APPROACH 

This  chapter  undertakes  the  derivation  of  a  negative-time  update  algorithm  built,  upon 
the  problem  formulation  of  the  previous  chapter.  This  approach  will  be  referred  to  as  a 
minimum  mean-square  error  (MMSE)  approach.  This  chapter  also  includes  a  simple  examp1  e 
to  display  the  use  of  this  approach,  and  compares  its  performance  with  other  methods  of 
incorporating  time-delayed  data. 


State  Evolution 


Referring  to  Figure  2-1,  the  illustrated 

T  = 

intervals  of  time  will  be  denoted  as 

tk  —  h~2 

(3.1) 

Ti  = 

tk-i  —  tk-  2 

(3.2) 

£3 

llt> 

tk  ~  tk-i 

(3.3) 

The  error  in  a  state  estimate  Xk\j  will  be  denoted  by  Xk\j  and  is  related  to  true  state  A*  by 

A*  =  Xk\j  +  Xky  (3.4) 

The  discrete-time  evolution  of  state  A*)*  6  TLn  over  the  time  interval  T  is  represented  as 

Afc  -  F(T)Xk-2  +  W(T)  (3.5) 


where  F  is  the  discrete  time  state  transition  matrix  and  W(T)  is  a  vector  representing  the 
effects  of  plant  noise  that  accrued  during  T.  Throughout  the  remainder  of  this  discussion,  it 
will  be  assumed  that  state  vec  tor  A*  =  [  xk  xk  )T  £c  one  (uncoupled)  dimension  in  a  two- 
or  three-dimensional  constant  velocity  filter.  The  state  transition  matrix  then  is 


F(T)  = 


1  T 
0  1 


(3.6) 


The  plant  noise  for  this  model  wil.  be  a  continuous-time,  white  noise  acceleration.  A  unit 
impulse  acceleration  S(t  ~  or),  with  ik-2  <  or  <  tk,  will  bring  about  the  response  in  state 
space  of  [(<*  —  <*)  1  ]3.  A  white  noise  acceleration  w(a)  where 


E{w{a)}  —  0,  E{w{a)w(f3)}  —  q  ■  8(a  —  (3) 


(3.7) 
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applied  throughout  the  time  interval  T  brings  about  (by  superposition)  the  plant  noise 

r** 


W(T)  =  I 

Jtk-2 


w(a)da 

The  time  projection  of  state  estimate  covariance  from  to  f*  is 


tk  ~  <* 
1 


(3.8) 


fl(W  =  F(T)Pk-W-iFT(T)+E  lj'\  [  ‘‘  j  “  w(a)iaf^  h  j 


1 T 


w(/3)dp\  (3.9) 

(3.10) 

(3.11) 


or 


with 


Pk\k- 2  =  F(T)Pk-2\k-iFT(T)  +  Q(T) 


Q(T)  = 


r  1 


|  T2  T 


where  q  is  a  scalar  parameter  in  units  of  (length2)/ (time3).  This  scalar  quantifies  the  un¬ 
certainty  in  the  (constant  velocity)  model  of  target  dynamics  and  will  be  referred  to  as  the 
“variance  of  acceleration  error.”  Thus,  when  the  time  interval  T  is  broken  into  two  segments, 
Ti  and  Ta,  the  state  transitions  are 


**-*  \  h-i  ~  « 


=  F(Tx)Xk-2  +  / 

Jtk 

Xk  =--  F(T2)Xk-i  +  £  [**" 


1 

tk-a 


w(a)da 


w(a)da 


and  the  projected  covariance  for  the  segmented  time  interval  is  given  by 

ft— i|k-2  =  ^(ft)ft-2|*-jFT(F,)  +  Q(ft) 
ft  |*-i  =  ftft)ft-m_,FT(ft)  +  Q(7y 

where 


Q(Tx)  = 
Q(T2)  = 


57? 

\T\ 

\T1 

Tx 

\n 

\n 

\n 

t2 

(3.12) 

(3.13) 

(3.14) 

(3.15) 

(3.16) 

(3.17) 


Filter  Update  At  ik  With  Zk 

When  the  measurement  Zk  arrives  at  the  filter  for  processing,  the  filter  does  not  know 
of  the  existence  of  the  delayed  measurement  Zk-j.  Therefore,  the  filter  will  proceed  with  a 
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standard  Kalman  filter  update  processing  using  T  =  tk  —  <*_3,  since  <*_3  was  the  time  tag  on 
the  last  measurement  that  it  received.  The  first  stage  of  processing  will  project  the  previous 
state  estimate  which  is  comprised  of  the  actual  state  less  the  estimation  error 


using  the  standard  time  projection  formula 

Xk\k-2  =  F(T)Xk-2\k-2 

Equation  (3.5)  shows  that  the  error  in  the  projected  state  estimate  it 

—  —  t*h  f  —  n  1 


**n-,  =  i».i + j‘‘^  r‘  j  “ 


u>(csr)d< 


(3.18) 

(3.19) 

(3.20) 


Standard  extended  Kalman  filter  processing  then  updates  the  state  estimate  using  gain  Kk 
and 

Xw  =  Xklk.2  +  Kk(Zk  -  HkXklk „3)  (3.21) 

Using  (3.20)  in  (3.21)  and  isolating  Xk  to  Xk\k  shows  that 


(I  -  KkHk)F(T)Xk-2\k-2  -  Kkvk 
+(/  -  KkHk)  /“"  f  **  7  “  1  w(a)da 

+(/  ~  KkHk)  ftk 

J*k-1 


tk-a 

1 


w(a)da 


(3.22) 

(3.23) 


To  allow  for  the  contingency  that  a  negative-time  update  may  be  necessary  at  the  time  of 
update  -X*_3|*_3  — »  it  will  be  convenient  to  define  the  quantity 


s*  t  E{xk[k  r 

Jtk-1 


tk-P 

1 


I  <v(,P)dP} 


(3.24) 


This  quantity  is  the  covariance  between  the  state  estimation  error  Xk\k  after  the  standard 
update  at  tk  and  the  plant  noise  error  that  accumulated  over  a  portion  (T2)  of  the  previous 
time  interval  T.  This  covariance  is  now  identifiable  by  substituting  (3.23)  into  (3.24).  The 
covariance  is  given  by 


Sk  =  £  j(/  -  KkHk)F(T)Xk„2\k.2  jf  ^ 


tk-P 

1 


w(/3)dp 


+(I  —  KkHk)  h  ^  w(a)daj *  tk  ^ 


tk-P 


+ (I~KkHk)jt^  tk  l°t  w(a)dajt ^  h  ~ 


w(p)d(3 


t v(P)dp 


tk-P 


-xkvk  r  h 

Jtk-i  1 


w(P)dp 


) 


(3.25) 
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Distributing  the  expected  value  through  (3.25)  gives 


Sk 


=  (I-KkHk)F(T)E{xk-2lh-2  l" 

1 


'*  t  |  »(w 


} 


+(/  -  KtH„)E  {£;  [y  w(a)daj'^  h  j  w(P)dfl\ 


+(/  -  |  £  [  '*  7  “  ]  j£,  [  **  1 


iT 


M/W 


(3.26) 


Eliminating  expectations  of  uncorrelated  random  processes  and  evaluating  the  single  remain¬ 
ing  term  yields 

Sk  =  (I-KkHk)Q(T2)  (3.27) 

which  is  the  sought-after  covariance  that  will  appear  in  subsequent  formulations. 


Negative-Time  Update  Algorithm 


When  the  occurrence  of  a  delayed  measurement  Zk-\  having  error  v*_i  ~  jV(0,  Rk-X) 
comes  to  light  after  the  state  estimate  update  Xk\k,  a  second  state  estimate  update 
using  an  as-yet  undetermined  Kalman  gain  Kk+  at  tk  will  be  made.  The  form  to  be  used  is 

Xk\k+  =  Xklk  +  Kk+{Zk_x  -  (3.28) 

Note  that  the  term  F~l(T2)Xk\k  is  used  to  compute  the  innovation  of  the  update  at  time 
tk„ i  while  keeping  the  time  of  the  estimate  Xk\k+  at  tk.  Using  Zk-\  ~  Hk„ \Xk-i  +  Vfc_i  and 
using  (3.13)  in  (3.28)  yields 


Xklk+  =  [I  -  Kk+H^F-'p^X^  +  Kk+Hk_xF-\T2)  f" 
The  modified  Kalman  gain  Kk+  is  now  found  by 


tk~a 

1 


w(a)da  -  Kk+vk_x 
(3.29) 


1.  Using  (3.29)  to  express  />*)*+  =  E{Xk\k+X^k+) 

Q 

2.  Setting  -r— ; — trace(Pk |*+)  to  zero 

oKk+ 

3.  Solving  for  Kk+ 
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In  so  doing,  at  is  useful  to  note  that  [4] 

4jtrace(BAC)  =  BTCfr  (3.30) 

OA 

•~trace(ABAT)  =  A(B  +  BT)  (3.31) 

OA 

The  covariance  P*|*+  as  £{^*1*+^*+}  is  found  from  (3.29)  to  be 

PM.+  =  -Ki+tft-, F-,(7i)]A|*|/ -  *»+#»-, F-,(7V)]T 
+[/  -  f,+ H„.IF-l(T1)]StF-T(T2)H?_lKl+ 
+K^Ht.lF-'(T1)Sf[I  -  Kki.Hi.lF-'(T,))T 
-t-KHiHi.lF-'(T1)Q(T,)F-T(T,)Ht1  +  «*-i]A'J+  (3.32) 

where  F~T  =  (F-1)T.  Setting  the  derivative  of  the  trace  to  zero  gives 

^~—trace(Pi,\k+)  =  -2ft|tF-T(r1)Pf_1 

+  2F*+[ff*-, F->(Ta){P*n  -  S,  -  Sf  +  Q(Ti))F-T(T1)Hl_l 
+  A_J  + 2S*F-r(7J)ffJ.1  =.  0  (3.33) 

Solving  for  the  modified  Kalman  gain  for  a  late  measurement  Kk+  yields 

Kk+  =  (P*|»  -  S„)F-T(T,)Hl,  x 

[ffl-,F-I(ra){P,|l  -s„-sl  +  <3(r1))F‘r(r1)pr_1  +  (3.34) 

Manipulating  (3.32)  into  the  form 

P„ |*+  =  F»+[ff4.,F-'(r1)(F.n  -  5*  -  Sf  +  <J(TJ)}F-T//r.,  + 

+F*I*  -  Kn.Ht-,F-l(T,)lPk |*  -  Sj] 

-(F*,»  - SJF-^TjHltKl  (3.35) 

and  substituting  (3.34)  into  the  first  Kk+  term  of  (3.35)  gives 

PkfA+  =  (/  -  F.+ff,-1F-‘(r1)]F,|1  +  KtkHt-,F-‘(T,)SZ  (3.36) 

which  together  with  (3.34)  and  (3.28)  completes  the  negative-time  update.  This  negative¬ 
time  update  algorithm  will  be  referred  to  as  the  MMSE  approach. 

Simple  Numerical  Example 

To  compare  the  MMSE  negative-time  approach  with  the  standard  Kalman  filter  formu¬ 
lation,  when  no  negative  times  are  encountered,  a  simple  numerical  example  is  presented 
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that  uses  the  Kalman  filter  for  one  coordinate  of  a  target  track.  The  filter  parameters  are 


Xk  = 

Fk(T)  = 

Hk  = 

R„  = 

q  = 


Xk 

1  T 
0  1 

1  0  J 

144.0 

5.0 


(3.37) 

(3.38) 

(3.39) 

(3.40) 

(3.41) 


where  the  plant  noise  is  computed  using  (3.11).  A  useful  parameter  for  characterizing  this 
tracking  filter  is  the  tracking  index  in  [3]  given  by 


He 


T3—  =  0.186 

Rk 


(3.42) 


The  tracking  index  is  often  useful  for  characterizing  tracking  filters  since  it  defines  an  a,  0 
filter  having  gains  equivalent  to  the  steady-state  Kalman  gains.  For  this  example  the  initial 
covariance  was  set  at 

"  79  32 
32  26 


P 3JO  = 


(3.43) 


A  baseline  case  is  now  tabulated  with  four  measurements,  Zk,  k  =  1,2, 3,4,  all  arriving  on 
time  at  1  sec  intervals.  Only  the  error  covariance  updates  are  tabulated.  The  equations  used 
are  the  standard  Kalman  filter  equations  given  by 


p*i*-i  =  P(T)pfc_I|fc_1pT(r)  +  Qk-i(T) 

(3.44) 

Kk  =  Pk\k-iHl[HkPk\k-iHt  +  Rk)'1 

(3.48) 

Pk\k  =  ( I-KHk)Pk\k-x 

(3.46) 

The  results  are  presented  in  Table  3-1 .  Table  3-2  shows  the  evolution  of  the  error  covariance 
matrix  for  the  case  of  measurement  Z3  missing  altogether.  The  final  measurement  Z4  retains 
the  subscript  4  to  conform  to  the  indexing  conventions  used  in  the  previous  section.  Note 
that  the  time  interval  T  increases  to  2  sec  and  that  the  elements  in  P4\2  are  quite  large.  The 
Kalman  gain  K4  is  also  high,  causing  the  trace  of  P4 14  to  increase  significantly  over  that  of 
the  baseline  case. 

The  negative-time  update  formulations  are  now  to  be  employed.  It  is  supposed  that  after 
the  update,  resulting  in  P4\4  with  Z3  missing,  has  taken  place  (but  before  any  subsequent 
updates),  Z3  arrives  late.  First,  the  covariance  Sk  is  found  by  (3.27)  to  be 


S4  =  (I  -  K4H4)Q(1)  = 
Using  P<|4  and  S4  in  (3.34)  yields 


0.6100  0.9000 
2.2000  4.5000 


0.3375 

0.0591 


(3.47) 


(3.48) 
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Table  3-1.  Baseline  Result  With  All  Measurements  On  Time 


k 

T 

Pk\k 

-  x 

Kh 

Pk\k 

1 

1 

'  163.67 
60.50 

60.50  ' 
31.00 

'  0.5320  ' 
0.1966 

76.60 

28.32 

28.32  ' 
19.10 

2 

1 

'  154.00 
49.92 

49.92  ' 
24.10 

‘  0.5168  ' 
0.1675 

74.42 

24.12 

24.12 

15.74 

3 

1 

‘  140.07 
42.36 

42.36  ‘ 
20.74 

0.4931 

0.1491 

71.00 

21.47 

21.47  ' 
14.42 

4 

1 

‘  130.04 
38.40 

38.40 

19.42 

0.4745  ' 
0.1401 

‘  68.33 
20.18 

20.18 

14.04 

Table  3-2.  Filter  Result  With  Measurement  Z3  Discarded  Altogether 
k  T _ Pk\k-i _ j£k _ Pk]k _ 


'  163.67  60.50  ' 
60.50  31.00 

'  0.5320  ' 
0.1966 

'  76.60  28.32  ' 
28.32  19.10 

‘  154.00  49.92  ' 
49.92  24.10 

'  0.5168 
0.1675 

'  74.42  24.12  ' 
24.12  15.74 

— 

— 

— 

'  247.20  65.60  ' 
65.60  25.74 

'  0.6319  ' 
0.1677 

'  91.00  24.15  ‘ 

24.15  14.74 

*  ■ 
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Table  3-3.  Traces  of  Covariances  P4 |4 


Case 

Trace 

All  measurements  on  time 

82.376 

Measurement  Z3  missing 

105.734 

Z3  processed  negative-time 

82.376 

and  using  this  result  in  (3.36)  yields 


^4|4+  = 


68.33  20.18 
20.18  14.04 


(3.49) 


Table  3-3  compares  the  traces  of  these  three  final  P4 14  matrices.  The  performance  of  the 
MMSE  approach  results  in  a  nearly  full  recovery  of  the  filter  covariance  achieved  when  the 
measurement  Z3  arrived  on  time. 


Negative-Time  Updates  and  Constant  Velocity  Tracking 


To  compare  the  MMSE  negative-time  update  with  other  methods  for  processing  mea¬ 
surements  with  a  negative-time  increment,  a  scenario  is  considered  in  which  there  are  two 
ships  transmitting  constant  velocity  track  information  periodically,  and  one  measurement  ar¬ 
rives  after  a  delay  that  exceeds  the  sample  period,  as  it  does  in  the  example  of  the  previous 
section.  This  situation  results  in  an  out-of-sequence  measurement  for  which  a  negative-time 
update  technique  must  be  utilized.  A  program  was  written  to  evaluate  the  effects  on  filter 
covariance  of  this  single  negative-time  update  when  tracking  a  constant  velocity  target  with 
various  negative-time  update  formulations.  Each  evaluation  includes  a  computation  of  filter 
covariance  similiar  to  the  numerical  example  of  the  previous  section. 

To  assess  the  performance  of  negativetime  update  techniques  for  various  tracking  sit¬ 
uations,  the  standard  deviation  of  the  measurement  errors  will  be  fixed  at  25  m,  while  the 
standard  deviation  of  the  acceleration  errors  for  the  filter  will  be  varied  from  2  to  25  m/s  t. 
Sample  periods  of  0.25,  0.5,  1.0,  and  2.0  sec  will  be  considered.  For  this  study,  each  of  two 
sensors  provides  data  at  a  sample  interval  of  2T,  with  the  data  from  one  sensor  skewed  by 
T  relative  to  the  other  sensor.  Thus,  the  measurements  of  the  target’s  position  are  recorded 
periodically  with  period  T  and  the  filter  will  achieve  steady-state  conditions.  The  filter  is 
assumed  to  be  in  steady-state  conditions  when  a  measurement  from  the  remote  sensor  is 
delayed  by  more  than  T,  and  a  negative-time  update  is  needed.  The  starting  covariance 
matrix  for  each  computation  is  derived  from  the  prevailing  tracking  index  pc  for  each  case. 
For  the  constant  velocity  tracking  filter  defined  by  (3.11)  and  (3.37)-(3.39),  the  steady-state 
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error  covariance  is  given  in  [3]  as 

Pss  =  Pk-l\h-\  ~  Rk 


where 

fi  =  6  -  3a  -  >/3a3  -  36a  +  36  (3.51) 

.  _  P  _T*1 

l-o  Rk 

With  Pss  of  (3.50)  as  the  starting  point,  the  usual  Kalman  filter  equations  were  used  to  up¬ 
date  the  covariance  matrix  using  a  sample  period  of  2T  to  represent  a  skipped  measurement. 
After  the  filter  update  with  sample  period  of  2T  is  completed,  the  delayed  measurement 
arrives  and  this  information  can  be  used  to  correct  the  state  estimate  and  associated  error 
covariance.  Five  techniques  for  processing  the  delayed  measurement  are  considered.  The 
first  three  are  those  illustrated  in  the  previous  section  with  the  simple  numerical  example. 
The  first  technique  (no  late  measurements)  is  an  ideal  case  in  which  the  measurement  is  not 
delayed  at  all.  The  second  technique  (discard  late  measurements)  is  put  forth  as  a  worst 
case  -  the  results  one  obtains  from  not  trying  to  cope  with  delayed  data  at  all.  The  third 
technique  (MMSE)  is  offered  as  a  practical  best  case  -  the  results  one  obtains  when  the  trou¬ 
ble  is  taken  to  fully  implement  an  optimized  negative-time  update  formulation.  The  next 
two  techniques  illustrate  the  benefits  and  limitations  of  ad  hoc  techniques  for  negative-time 
updates.  The  fourth  technique  includes  plant  noise  Qk  during  the  time  update  to  the  late 
measurement  and  during  the  time  update  back  to  the  current  time.  The  fifth  technique  sets 
the  plant  noise  Qk  to  zero  during  the  time  update  to  the  late  measurement  and  during  the 
time  update  back  to  the  current  time. 

The  five  techniques  are  compared  using  a  plot  of  the  normalized  trace  of  the  error  co- 
variance  given  by 

flroce^j  (353) 

versus  the  standard  deviation  of  the  acceleration  error.  The  results  are  given  for  each  time 
period  in  Figures  3-1  through  3-4.  Using  all  late  measurements  discarded  as  a  worst  case 
senario,  it  is  seen  that  the  technique  using  the  plant  noise  throughout  the  negative  update 
breaks  down  for  relatively  low  levels  of  acceleration  error.  Breakdown  of  a  negative-time 
update  algorithm  will  be  used  to  denote  the  condition  at  which  the  error  covariance  at 
time  tk  is  increased  by  processing  a  time-delayed  measurement  at  time  —  A.  However, 
the  technique  that  has  the  plant  noise  set  to  zero  during  the  negative  update  seems  to  be 
sufficient  for  shorter  time  periods,  but  becomes  unrealistic  for  longer  time  periods  in  that  the 
trace  of  the  error  covariance  is  lower  than  the  optimal  case  where  all  measurements  arrive 
sequentially  (no  late  measurements).  It  is  also  noticed  that  the  performance  of  the  MMSE 
negative  update  closely  matches  that  of  the  sequential  case. 


(3.52) 


-  | 

£  -  0 

T  2(1  -a)T* 


(3.50) 
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STANDARD  DEVIATION  OF  ACCELERATION  ERROR 

Figure  3-1.  Negative-Time  Update  Technique  Performance 


STANDARD  deviation  of  acceleration  error 


Figure  3-2.  Negative-Time  Update  Technique  Performance 
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STANOARO  DEVIATION  OF  ACCELERATION  ERROR 

Figure  3-3.  Negative-Time  Update  Technique  Performance 


STANDARD  DEVIATION  OF  ACCELERATION  ERROR 


Figure  3-4.  Negative-Time  Update  Technique  Performance 
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CHAPTER  4 
CEC  APPROACH 

This  chapter  interprets  the  ad  hoc  CEC-extended  Kalman  filter  updating  technique  as 
specified  in  [1].  Analyse  similar  to  that  in  the  previous  chapter  is  included  for  examination 
of  the  technique’s  performance. 

CEC  Specification  For  Negative-Time  Updates 

Assuming  the  same  situation  as  in  the  previous  chapter,  a  measurement  arrives  late  with 
time  tag,  f*_i,  just  after  the  extended  Kalman  filter  update  has  been  made  for  a  measurement 
with  time  tag,  tk ■  This  results  in  the  time  interval 

T  =  4-i  -  h  (4.1) 

From  this  point,  noting  T  <  0,  the  standard  filtering  equations  are  used  to  make  a  back¬ 
wards  update  with  the  exception  that  the  plant  noise  is  set  to  zero,  (i.e.,  Q(T)  =  0).  The 
equations  for  this  segment  of  the  interpreted  CEC  negative-time  update  are  given  in  the 
following  equations. 

Time  Update: 


Xk-w-  =  F{T)Xk  i*-  (4.2) 

=  F{T)P^-Ft(T)  (4.3) 


Measurement  Update: 

+  Kk-i[ZK_i  —  h(Xk-.i\k-)]  (4.4) 

Pk- 1|*  =  [I  -  (4.5) 

Kt-1  =  Pk-nk-HLAHt-tPk-^-Hh  +  (4.6) 


After  the  state  has  been  moved  backwards  and  the  state  estimate,  has  been  calcu¬ 

lated,  the  state  estimate  is  moved  forward  with  the  plant  noise  reintroduced  by  using  the 
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absolute  value  of  the  negative-time  interval  to  project  forward  to  the  current  track  valid 
time,  ffc,  to  form  -X*)*.  This  segment  gives  the  time  updates  as 

X*|W  =  Fm)Xi-w  (4.7) 

ft  l*+  =  FdTDn-.l^dT))  +  Q(|T|)  (4.8) 

noting  that  now  <?(|r|)  is  nonzero  and  the  time  interval  |T|  >  0.  Having  predicted  for¬ 
ward  to  the  existing  track  valid  time,  tk ,  the  filter  resumes  normal  operations  on  the  next 
measurement. 


Simple  Numerical  Example 


A  simple  example  that  computes  the  covariance  matrix  is  presented  for  comparison  with 
the  results  of  the  previous  example  in  Chapter  3.  The  following  parameters  were  chosen  to 
remain  consistent  with  the  previous  example.  Those  parameters  are  given  by 


X*  - 

X 

X 

k 

(4.9) 

Fk(T)  = 

'  1 

0 

T  ' 
1 

(4.10) 

Hk  =  | 

’  1 

°] 

(4.11) 

Rk  =  144.0 

(4.12) 

q  =  5.0 

(4-13) 

The  initial  covariance  matrix  is  set  to 


Polo 


79  32 
32  26 


(4.14) 


The  CEC  approach  is  used  to  calculate  the  error  covariance  matrix  using  (4.2)-(4.8)  as 
discussed  earlier.  The  final  error  covariance  matrix  results  are 


P||4+  - 


70.48  23.56 
23.56  19.30 


(4.15) 


Recalling  the  results  from  the  previous  example  found  in  Tables  3-1  and  3-2,  a  comparison  of 
all  four  cases  is  presented  in  Table  4-1.  It  is  noted  that  although  the  CEC  ad  hoc  technique 
made  a  significant  recovery  from  the  negative-time  interval,  a  nearly  full  recovery  was  not 
experienced  as  with  the  MMSE  approach. 
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Table  4-1.  Traces  of  Covariances  P4 j4 


Case 

Trace 

All  measurements  on  time 

82.376 

Measurement  Z3  missing 

105.734 

Z3  processed  negative-  time(MMSE) 

82.376 

Z3  processed  negative-time(CEC) 

89.777 

Negative-Time  Updates  and  Constant  Velocity  TVacking 


The  program  method  used  in  the  last  section  of  Chapter  3  was  used  to  evlauate  the 
performance  of  the  interpreted  ad  hoc  CEC  technique.  The  worst  case  scenario,  where  all 
late  measurements  are  discarded  and  the  sequential  case  where  all  measurements  arrive  on 
time,  were  compared  with  the  CEC  technique.  The  performance  was  again  evaluated  using 
a  graph  of  the  the  normalized  trace  of  the  filter  error  covariance  matrix  versus  the  standard 
deviation  of  the  acceleration  error.  Time  periods  of  0.25,  0.5,  1.0,  and  2.0  sec  were  used  in 
the  evaluations.  Figures  4-1  through  4-4  express  the  performance  of  the  three  techniques. 
The  ad  hoc  CEC  technique  does  not  always  perform  better  than  the  worst  case  scenario  of 
discarding  all  late  measurements.  For  longer  time  periods,  the  CEC  technique  was  found  to 
perform  poorly,  even  with  relatively  small  standard  deviations  of  acceleration  error,  while 
the  MMSE  technique  remained  stable  as  shown  in  the  previous  chapter. 

For  each  of  the  time  periods  used,  the  standard  deviation  of  acceleration  error  at  which 
the  CEC  technique  breaks  down  was  evaluated.  Using  (3.42),  the  corresponding  tracking 
indexes  were  found  in  order  to  characterize  the  break  down  regions  for  the  CEC  technique. 
Figure  4-5  shows  a  plot  of  the  tracking  indexes  at  which  the  CEC  technique  fails  versus  the 
sample  period.  For  a  given  sample  period,  the  CEC  technique  will  break-down  for  tracking 
index  values  greater  than  the  value  denoted  by  the  line. 
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- 1 - f  . . 1 - r 

1  —  NO  LATE  MEASUREMENTS 

Z  —  DISCARD  LATE  MEASUREMENTS 
3  —  CEC  TECHNIQUE 


TIME  PERIOD  —  r>.  25  SEC 


i _  _ _ -i _ i _ i _ 

O  S  1  O  IS  20 

STANDARD  DEVIATION  OF  ACCELERATION  ERROR 

Figure  4-1.  Negative-Time  Update  Technique  Performance 


o  ■ - 1 - 1 - 1 - 1 — 
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STANDARD  DEVIATION  OF  ACCELERATION  ERROR 

Figure  4-3.  Negative-Time  Update  Technique  Performance 


Figure  4 . 4.  Negative-Time  Update  Technique  Performance 
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Figure  4-5.  Tracking  Index  Break-Down  Region  for  the  CEC  Technique 
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CHAPTER  5 

SIMULATION  RESULTS 


The  objective  was  to  simulate  an  ownship  sensor  and  a  remote  ship  sensor  tracking  a 
common  manuevering  target  with  the  remote  ship  transmitting  its  tracking  data  to  ownship. 
This  simulation  enabled  the  performance  of  the  interpreted  CEC  negative-time  approach  to 
be  compared  with  the  derived  MMSE  approach  in  a  practical  situation  where  data  being 
transmitted  from  remote  ship  to  ownship  are  delayed. 


Simulation  Setup 

Figure  5-1  displays  the  scenario  of  the  ships  and  target  locations  for  this  simulation. 
The  ships  are  separated  by  15  km  along  the  xaxis  making  the  remote  ship  coordinates  with 
respect  to  the  ownship  reference  frame, (x,p,  z)  =  (15,0,0)  km.  To  simplify  the  simulation, 
motion  between  the  two  sensors  was  set  to  zero  and  the  earth’s  shape  was  modeled  to  be 
locally  flat.  The  flow  chart  in  Figure  5-2  is  a  basic  outline  of  the  program  operations. 
A  single  target  trajectory,  created  using  a  trajectory  generator  program,  was  used  for  the 
simulation.  A  profile  of  the  trajectory  is  found  in  Figure  5-3.  The  target  has  a  constant 
altitude  of  5  km.  The  target  travels  at  Mach  2  for  44  sec  making  two  6g  maneuvers  between 
12  and  21  sec  and  between  24  and  33  sec,  respectively.  When  the  target  is  not  maneuvering, 
it  maintains  a  constant  velocity.  The  trajectory  generated  is  in  the  form  of  a  data  matrix 
with  even  columns  representing  measurements  from  the  ownship  sensor  and  odd  columns 
representing  the  measurements  obtained  by  the  remote  ship  sensor.  The  odd  columns  are 
offset  by  (15,0,0)  km,  since  every  other  measurement  received  by  the  ownship  is  a  remote 
sensor  measurement.  Measurements  from  each  sensor  are  received  every  2  sec  resulting  in  a 
received  measurement  every  1  sec  when  the  data  from  each  sensor  are  skewed  in  time. 

The  same  trajectory  was  ULed  in  the  simulation  to  study  and  test  each  of  four  cases.  For 
the  first  case,  the  ownship  uses  its  own  measurements  plus  measurements  from  the  remote 
ship  to  track  the  target,  with  all  measurements  arriving  on  time.  The  second  case  is  similar  to 
the  first  case  except  all  the  measurements  received  by  the  ownship  from  the  remote  ship  are 
delayed  and  are  discarded.  The  third  case  resembles  the  second,  but  the  late  measurements 
from  the  remote  ship  sensor  are  processed  using  the  MMSE  negative- time  update  approach 
derived  in  Chapter  3.  In  the  lourth  case,  the  delayed  measurements  received  by  the  ownship 
from  the  remote  ship  are  processed  using  the  interpreted  CEC  approach  of  Chapter  4. 
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Figure  5-1.  Profile  of  Ship  and  Target  Locations 
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Figure  5-2.  Basic  Outline  of  Simulation 
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Figure  5-3.  Profile  of  Target  Trajectory 
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Result* 

The  simulation  output  for  each  track  filter  is  an  average  of  100  Monte  Carlo  experiments. 
The  Root- Mean- Square- Errors  (RMSE)  in  position  and  velocity  are  used  to  measure  the 
performance  of  the  negative-time  update  techniques  implemented  in  the  simulations.  In 
the  four  cases,  the  combined  measurements  received  from  the  two  multiplatform  sensors 
had  a  data  rate  of  1  Hz  with  a  standard  deviation  in  range  of  <rT  —  12  m,  and  standard 
deviation  of  angles  o*  =  <re  =  2.5  mrad.  The  arrival  of  measurements  from  the  remote 
sensor  are  delayed  slightly  more  than  1  sec.  The  results  of  the  first  three  cases  are  given  in 
Figure  5-4.  The  top  and  bottom  graphs  are  an  average  RMSE  of  the  position  and  velocity, 
respectively.  Comparing  the  first  two  cases  displays  the  improvement  experienced  with  the 
concept  of  CEC,  when  all  measurements  from  the  remote  sensor  arrive  on  time,  as  opposed 
to  using  only  the  ownship  sensor  measurements.  In  the  third  case,  the  late  measurements  are 
processed  using  the  MMSE  negative-time  update  approach  derived  in  Chapter  3.  It  is  seen 
that  the  result  is  almost  a  complete  recovery  from  the  late  measurements.  The  jagged  lines 
that  occur  during  the  maneuvers  for  the  case  where  all  measurements  arrive  sequentially  are 
the  result  of  receiving  data  from  multiple  sensors  with  varying  degrees  of  accuracy.  Since 
the  remote  sensor  is  closer  to  the  target  than  the  ownship  sensor,  as  shown  in  Figure  5-1, 
measurements  from  the  remote  sensor  are  better  than  the  measurements  from  the  ownship 
sensor.  Since  the  measurements  between  ownship  and  remote  are  skewed  in  time,  every  other 
measurement  being  utilized  will  result  in  a  better  estimate.  During  the  constant  velocity 
portions  of  the  track,  this  effect  goes  unnoticed,  but  during  maneuvers  where  the  mismatch 
of  motion  model  is  greater,  the  effect  is  seen.  Since  the  MMSE  approach  does  not  compute 
a  state  estimate  that  corresponds  to  the  measurement  time  of  the  remote  sensor,  RMSE 
were  not  computed  at  those  times,  and  hence,  the  lines  denoting  the  RMSE  for  the  MMSE 
approach  are  not  jagged.  The  results  comparing  the  RMSE  of  position  and  velocity  for  cases 
three  and  four  are  found  in  Figure  5-5.  Comparing  the  two  techniques,  it  is  seen  that  the 
CEC  approach  does  result  in  a  recovery  from  the  late  measurements,  but  the  results  are  not 
as  precise  as  the  recovery  experienced  from  the  MMSE  approach. 
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Figure  5-5.  Tracking  Results  for  Cases  3  and  4 
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CHAPTER  0 
CONCLUSIONS 


An  MMSE  approach  for  processing  time  delayed  data  was  developed  and  compared  with 
other  techniques  for  processing  time  delayed  data.  One  of  those  techniques  for  processing 
negative-time  data  uses  the  plant  noise  Qk  throughout  the  entire  negative-time  update  in¬ 
terval  (Q-on).  Analysis  of  the  Q-on  technique  showed  that  it  is  unstable  and  performs  very 
poorly.  Another  technique  involved  setting  Qk  to  zero  throughout  the  entire  negative-time 
update  interval  (Q-off).  Analysis  of  the  Q-off  technique  showed  that  it  performs  better  than 
the  Q-on  technique,  but  for  larger  levels  of  plant  noise  and  longer  time  intervals  the  filter 
covariance  is  “too  good”  in  that  it  is  smaller  than  that  achieved  by  sequential  processing 
of  the  data  (all  measurements  on  time).  The  interpreted  CEC  approach  was  found  to  a  be 
more  desirable  approach  to  negative-time  updates  than  either  the  Q-on  or  Q-off  technique. 
However,  when  the  MMSE  approach  was  compared  with  the  CEC  approach,  the  CEC  ap¬ 
proach  was  found  to  be  the  inferior  technique.  The  CEC  approach  may  have  problems  when 
used  with  filters  designed  for  highly  maneuvering  targets  at  low  data  rates,  while  the  MMSE 
approach  is  rather  stable  with  respect  to  levels  of  plant  noise  and  sample  periods. 

The  possibility  of  experiencing  delayed  data  causing  negative-time  intervals  has  added 
complexity  to  the  CEC  concept.  More  analysis  of  the  issues  associated  with  negative-time 
updates  is  desirable.  Issues  revolving  around  the  negative-time  update  that  need  to  be  stud¬ 
ied  include  status  of  target  maneuvers,  data  association  in  the  presence  of  clutter  and  false 
alarms,  and  an  on-line  technique  to  assess  the  value  added  if  a  negative-time  measurement  is 
utilized.  This  study  involved  only  a  constant  velocity  Kalman  filter.  Other  problems  could 
arise  if  other  tracking  techniques  such  as  multiple  model  filtering  are  used.  For  example,  at 
the  beginning  of  a  target  maneuver,  a  multiple  model  filter  may  undergo  a  model  change  just 
prior  to  receiving  negative-time  data.  In  this  situation  it  would  be  useful  to  have  a  method 
for  deciding  which  model  would  be  used  for  processing  the  negative-time  data.  Also,  an 
on-line  technique  is  needed  for  evaluating  the  need  and  benefit  of  using  delayed  data  during 
a  tracking  process.  In  a  situation  where  good  data  on  a  target  has  been  attained  recently, 
the  delayed  data  can  be  ignored,  while  on  the  other  hand,  if  very  little  data  on  the  target 
have  been  attained  recently,  the  delayed  data  may  be  critical  to  maintaining  a  good  track. 
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APPENDIX  A 

NOTE  ON  R  R.  CASTELLA’S  MULTI-SITE,  MULTI-SENSOR  TRACKING  FILTER 


A-l 


NSWCDD/TR-93/351 


The  proposed  tracking  algorithm  described  in  this  appendix  first  appears  in  the  Applied 
Physics  Laboratory  internal  memorandum  F2A-92-1-002  dated  16  January  1992.  Its  principle 
features  are  as  follows: 

(1)  A  discrete  time  Kalman  filter  with  state  vector 

X=  [  x  x  y  y  z  z]T 
and  a  constant  velocity  model  $(T)  is  used. 

(2)  Measurements  from  various  sites,  including  ownship,  are  received  for  processing.  A  given 
measurement  may  come  from  one  of  several  types  of  sensors.  One  may  be  from  a  radar  that 
can  measure  elevation  but  not  doppler,  whereas  another  may  come  from  a  radar  that  can 
measure  doppler,  but  not  elevation. 

(3)  Measurements  Mk  arrive  as  “raw  data”;  that  is,  individual  measured  parameters  such  as 
bearing  or  slant  range  each  appear  accompanied  by  a  time  stamp,  a  standard  deviation,  sen¬ 
sor  location  (as  latitude-longitude),  and  sensor  ship’s  heading,  pitch  and  roll  as  required.  The 
Kalman  updates  are  accomplished  one  measured  parameter  at  a  time  with  an  appropriate 
output  matrix  Hk  computed  for  each  measured  parameter  at  each  update. 

(4)  Each  Kalman  filter  within  the  Cooperative  Engagement  System  operates  from  the  point 
of  view  of  a  local,  stabilized  (xj,yi,  Zy)  coordinate  system.  Data  from  other  sensors  are  ref¬ 
erenced  to  a  remote,  stabilized  coordinate  system  ( x2 ,  y2,  z2)  with  a  known,  remote  latitude- 
longitude. 

(5)  The  Kalman  filter  with  its  state  X  G  7Z6  has  the  usual  covariance  matrix  P  6  H6*6 
representing  accumulated  measurement  errors  and  plant  noise.  There  also  is  a  measurement- 
noise-only  covariance  matrix  Ps  €  7l6*6  that  is  carried  to  facilitate  decision  processes  outside 
the  Kalman  filter  proper.  The  matrix  Ps  does  not  otherwise  affect  the  operation  of  the  filter. 

The  Kalman  filter  equations  used  are  stated  below.  The  notation  used  here  is  X*|*_  j 
for  projected  state  (as  opposed  to  6£.  used  originally)  and  X*|*  for  filtered  state  (as  opposed 
to  6£.  in  the  original).  1 


=  *(r)x*_! 

(A.l) 

=  *(T)Pk^k-MT)T  +  Qk 

(A.2) 

=  $(T)Ps,k-i\k-i$(T)T 

(A.3) 

Kk 

=  Pk\k^Hl{HkPk\k-XHT +  Rk)~x 

(A. 4) 

=  X*,*-!  Kk[Mk  -  A/p(X*|*_i)] 

(A.5) 

Pk\k 

=  (I-KkHk)Pk 

(A.6) 

lrThe  original  document  treats  state  as  £(0),  the  state  at  track  initiation,  plus  a  summation  of  incremental 
corrections  Sj^.  This  convention  is  not  used  here, 
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Ps,k\k  =  (/  -  KkHk)Ps,k\k-i(I  -  KkHk)T  +  KkRkKj  (A.7) 

The  specific  means  by  which  these  Kalman  filter  equations  are  employed  to  effect  track 
filtering  will  now  be  explored. 


Plant  Model  and  Plant  Noise 


Since  the  Cooperative  Engagement  Processor  (CEP)  tracking  algorithm  is  a  constant- 
velocity  Kalman  filter,  the  transition  model  for  the  chosen  state  vector  X  is  given  by 


$(T)  = 


1  T  0  0  0  0  ' 

0  1  0  0  0  0 

0  0  1  T  0  0 

0  0  0  1  0  0 

0  0  0  0  1  T 

0  0  0  0  0  1. 


(A.8) 


where  T  is  the  time  increment  between  updates. 


In  this  particular  realization,  measurements  Mk  are  each  scalars  that  reflect  some  pre¬ 
dicted  function  of  the  state  vector  MP(X),  rather  than  just  one  of  the  elements  of  X.  The 
‘innovation’  (the  bracketed  term  of  (A.5)  above)  is  now  the  difference  between  measurement 
Mk  and  ‘predicted’  Mp(Xjt|*_i).  Thus,  the  measurement  matrix  Hi  becomes  a  row  vector  of 
partial  derivatives  given  by 


evaluated  at  X^ijt—x . 


dMp  dMp 
dx  dy 


dMp 

dy 


dMp 

dz 


dMp  ' 
dz 


(A.9) 


The  plant  noise  Qk  is  represented  by  the  formulation 

Qk 


q«  = 


'2x2  — 


0 

0 


r 

qi 

02x2 

02X2 

= 

02x2 

q2 

02x2 

.  02x2 

02X2 

qa  . 

r  t 3 

ji2  - 

=  9. 

_ i 

Y 

T 

,  i  = 

i  =  1,2,3 


(A.10) 

(All) 

(A.12) 


Each  scalar  9,  is  a  plant  noise  spectral  density  in  units  of  length2/ time3  =  acceleration2/Hz. 
This  parameter  will  be  supplied  based  upon  a  separate  determination  of  the  existence  of 
a  target  maneuver.  Further  details  of  maneuver  detection  and  acceleration  estimation  are 
not  given  in  the  memorandum  under  discussion  and  must  be  treated  in  further  commentary 
notes. 
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Coordinate  Systems 


With  the  role  of  the  CEP  tracking  algorithm  being  to  accept  measurements  from  many 
different  sensors  at  many  different  locations,  coordinate  conversion  plays  a  very  important 
part  in  its  formulation.  This  section  will  set  out  coordinate  system  notation  conventions, 
drawings  of  the  various  coordinate  systems,  and  the  relationships  used  to  express  measure¬ 
ments  that  are  made  in  one  system  in  terms  of  another. 


Figure  A-l  illustrates  worldwide  coordinate  systems  expressed  in  longitude  L  (positive  L 
=$■  East  from  Prime  Meridian),  latitude  p  (positive  p  =>  North  from  the  Equator),  and  in  the 
earth-centered  rectangular  coordinate  system  (xo,  yo,  Zo)-  Any  vector  Ri  can  be  expressed  as 
an  earth-centered  vector  (xo,yo,z0)  in  terms  of  Li,  p,',  and  its  length  |/2,|  by  the  well-known 
relationship 


Xo 

ccs  Li  cos  p, 

= 

yo 

sin  Li  cos  pi 

.  Z0  . 

sin  pi 

(A.13) 


In  Figure  A-l,  the  head  of  vector  Ri  is  shown  as  the  origin  for  a  ‘stabilized’  coordinate  system 
(x, that  is,  one  for  which  x  is  local  East,  y,  is  local  North,  and  £<  is  local  vertical.  In 
the  CEP  tracker,  R\  points  to  ownship  and  state  estimate  X  is  maintained  in  the  (xi,  y\,  Zi) 
coordinate  system.  Measurements  received  from  another  ship  located  at  (/,2,p2)  and  at  the 
head  of  vector  R2  may  have  been  made  in  either  stabilized  local  coordinates  as  (x2,  y2,  z 2)  or 
in  ‘deck  coordinates’  as  (x2",  y"' ,  z’£)  having  the  same  origin. 


Figure  A-2  illustrates  the  relationships  among  earth-centered  coordinates  (x0,  yo,  Zo), 
local  stabilized  coordinates  (xj,yi,zi),  and  remote  stabilized  coordinates  (x2,y2,  z2).  An 
ownship  sensor  makes  measurements  on  a  target  in  terms  of  vector  Tx,  while  a  remote  sensor 
makes  measurements  on  the  same  target  in  terms  of  vector  T2-  The  vector  displacement 
between  remote  and  local  coordinates  in  earth-centered  (x0,yoi  ^o)  is 


zkffo  —  Ri  ~  Ri 

The  well-known  matrix 


\Ri\  cos  fi\  cos  L\  —  j/?2J  cos  ^2  cos  L2 
|/?i  |  cos  p  1  sin  L\  —  |/?2|  cos  p2  sin  L2 
|Ri|  sinpi  —  |A2|  sin  p, 


(A. 14) 


—  sin  Li  cos  L\  0 

—  sin  p,  cos  Li  sin  p,  sin  L,  cos  p, 
cos  p,  cos  Li  cos  pi  sin  T,  sin  p, 


(A.15) 


rotates  any  vector  (xo>2/o>  *o)  into  the  i-th  coordinate  system  from  the  earth- centered  system. 
Coordinate  rotation  in  the  opposite  direction  is  achieved  via  the  inverse  matrix 

=-  [M(,0)]-1  =  (A. 16) 


which  is  equal  to  the  transpose  since  these  matrices  are  unitary.  Thus,  the  displacement 
between  local  (xi,t/i,Zi)  and  remote  (x2,y2,  z2)  in  local  coordinates  is 

A  Rx  =  M(,0)A/?0  (A. 17) 
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Figure  A-l.  Longitude  L,  Latutude  /i,  and  Earth- Centered  Unit  Vectors  (x0,  yo,  Zq) 


Figure  A-2.  Local  and  Remote  Coordinate  Systems  of  the  CEP  Tracker 
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Figure  A-3.  Converting  Stabilized  Coordinates  to  Deck  Coodinates 

Hence,  any  attribute  of  a  target  measured  in  terms  of  vector  T2  by  a  remote  sensor  at  (L2,  H2) 
as  stabilized  {x2,y2,z2)  is  expressable  in  terms  of  the  ownsliip  system  state  vectors 


x2 

1 

'  Xi  ' 

2/2 

=  A/<20>M<01> 

yi 

*2  . 

V 

.  zi . 

H  Aili 


(A. 18) 


Some  ship  sensors,  such  as  the  SPY-1,  wiil  produce  their  measurements  in  stabilized 
coordinates  ( x  as  East,  etc.),  while  others,  such  as  the  SPS-49,  report  in  ‘deck  coordinates’ 
(ym  is  toward  the  bow  of  the  ship).  Figure  A-3  illustrates  the  three  angles  by  which  a 
stabilized  system  is  rotated  into  a  deck  system.  The  Navy  convention  for  coordinate  rotation 
differs  from  the  mathematical  convention  in  the  order  and  in  the  sign  convention  of  the 
individual  rotations.  In  the  Navy  convention,  positive  angles  represent  clockwise  rotations. 
In  going  from  the  stabilized  system  (x„y,-,z,)  to  the  deck  oriented  system  (x'/'.  y"1,  z"'),  one 
first  rotates  about  z,  by  heading  hi  (which  is  positive  for  clockwise  from  North)  into  (x$,  y-,  zQ 
by  way  of 


(A.19) 


Secondly,  one  rotates  about  x(  by  the  pitch  angle  pi  (which  is  positive  for  the  bow  going 
down)  to  get  into  (xf,  y",  z")  by  way  of 


■*51 

cos  hi 

—  sin  hi  0  " 

'  Xi  ' 

t/i 

-  H{i)Xi  = 

sin  hi 

cos  hi  0 

Vi 

Z 

*1  J 

0 

0  1  _ 

.  Zi  . 

r  *" 1 

y'i 

=  P^X'i  = 

L  z"  J 

1  0  0 
0  cos  pi  —  sinp, 
0  sin  pi  cos  pi 


'  x'i  ' 

y[ 

.  zl . 

(A. 20) 


A-7 


NSWCDD/TR-93/351 


Finally,  one  rotates  about  y"  by  roll  angle  r i  (which  is  positive  for  roll  to  portside)  into  the 
final  deck  system  (i-*,  y(",  z"')  by  way  of 


cos  r, 

0 

sinr< 

r  *r  1 

Vi 

s tl' 

=  R(,)X"  = 

0 

1 

0 

1  *r  - 

—  sinr,- 

0 

COS  Tj 

z! 

L 

(A.21) 


Thus  measurements  made  by  an  ownship  sensor  in  deck  coordinates  are  expressed  in 
ownship  stable  coordinates  as 


r"'  1 
X1 

'  Xi  ' 

y'x 

= RWpWHW 

Vi 

JH 

*-1  J 

.  rl  . 

(A.22) 


and  a  remote  sensor  using  deck  coordinates  makes  measurements  expressed  in  ownship  stable 
coordinates  as 


t’" 

Ml 

Vi 

Lb  J 


rWpWhMmWmI 01) 


*1 

Vi 

Zl 


-f  Ai?i 


(A.23) 


For  convenience, 


the  following  products  of  coordinate  rotation  matrices  will  be  designated. 


U  =  RWPWH(1) 

(A.24) 

N  =  M<20>M<01> 

(A.26) 

V  £  RWpWHWN 

(A.26) 

Measurement  Functions  Mp  and  Matrices  H 

The  CEP  tracking  algorithm  is  configured  to  handle  four  kinds  of  measurements  from 
two  different  kinds  of  sensor  (radars),  neither  of  which  alone  can  make  all  four  kinds  of  mea¬ 
surements.  One  kind  of  radar,  the  SPY-1  type,  produces  measurements  in  local  stabilized 
coordinates  whereas  other  types,  bearing-only  sets,  produce  their  results  in  deck  coordi¬ 
nates.  Moreover,  a  given  processor  my  be  resident  on  a  ship  having  either  type  of  set.  The 
measurement  types  are 

Bearing 

Slant  range 

Elevation 

Doppler 
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In  order  to  process  the  parameter  M*  characteristic  of  each  kind  of  measurement,  the  fol¬ 
lowing  steps  must  be  undertaken. 


(1)  Evaluate  Mp(X*|*_i)  for  predicted  parameter. 

(2)  Evaluate  the  partial  derivatives  needed  to  complete  matrix  /f£. 

(3)  Execute  the  Kalman  filter  formulas  (A-l)  through  (A-7). 


The  formulas  needed  to  carry  out  Steps  (1)  and  (2)  will  now  be  discussed. 


Formulas  for  MP(X.)  and  c  =  x,x,y,y,z,z,  for  each  measurement  type  are  given 


dM0 


in  Tables  A-l  through  A-4.  Table  A-l  gives  MP(X)  and  — for  processing  local  stabi¬ 
lized  coordinates  (a:lv y\,zx),  while  Table  A-2  gives  those  results  for  local  deck  coordinates 

3  M 

(*T,nr.*r).  Table  A-3  gives  MP(X)  and  for  processing  remote  stabilized  coordinates 

(xi,y2,  Z2),  while  Table  A-4  gives  those  results  for  remote  deck  coordinates  (x™,  j4",  z'f).  As 
an  example  of  these  formulas’  derivations,  the  bearing  measurement  will  be  treated  in  detail. 


Bearing  in  Local  Stabilized  Coordinates 


While  bearing  may  be  defined  in  terms  of  the  tangent  function  tan  Bx 


the  inverse  B\ 


=  tan-  (S) 


(xA  ■ 

—  I ,  in  usi 

\yij 


using 


the  programmer  must,  of  course,  make  provisions  for  cases 


where  yx  <  0.  Since  these  provisions  involve  adding  a  constant  to  Bx(xx,yx),  the  formulas 
for  partial  derivatives  need  no  such  modifications. 


Using  the  general  formulation  for  the  derivative  of  tan 


-1 


du  dv 

vdTx~ute~x 

u2  +  v2 


With  u  =  z  1  and  v  =  yx  it  follows  immediately  that 

JL  r  Mi 

dxx  1  x2  +  y2 


and  that 

—B  ~Xl 
dyx  1  x2  +  y2 


(A.27) 


(A.28) 


(A.29) 
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Bearing  in  Local  Deck  Coordinates 
Now  bearing  is  defined  by 

tan5i"=37  (A.30) 

Vi 

Using  the  coordinate  rotation  matrix  U  =  [£/,j]  of  (A.24),  x'"  and  yj"  are  found  to  be 


if  =  [l  0  o]t/ 


*1 

Vi 

z\ 


—  x\U\\  +  y\U\i  +  z\U\% 


and 


y'"  —  [  0  1  0  ]  t/ 


yi 

L  J 

=  +  J/lC^M  +  Z1U23 

Again  using  (A. 27)  with  u  =  x'"  and  v  =  y“'  the  derivatives  become 

_  y?Un  -  *?Un 

dxi  (x?y  +  WY 

with  a°d  following  similarly. 


Bearing  in  Remote  Stabilized  Coordinates 


In  this  instance 


tan  82  —  — 

y2 


where  a?2  and  y2  are  found  from  (A. 18)  and  (A.25)  as 

[l  0  0  ]  N 


x2  = 


y2  = 


ii  ■+■  A  Rzi 
yi  +  A/lyl 
L  z\  +  A  RA 

•Nn(^i  +  A/1,1)  +  Nu(yi  +  A/lyj)  +  N\3(zi  +  A/2,i) 

zi  +  A/lxi 
yi  +  A/?yi 

L  *1  +  A/i,i  J 

■^21(^1  +  A/?,i)  +  •N32(yi  +  A Ry2)  -f  #23(21  -f  A/?,i) 


[  0  1  0  ]  N 


(A.31) 

(A.32) 

(A.33) 

(A.34) 

(A.35) 


(A. 36) 

(A.37) 

(A.38) 

(A.39) 

(A.40) 
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Derivatives  also  follow  from  (A.27)  with  u  =  x2  and  v  =  ya  so  that 

dBj _ yjJVn  —  ijiVji 

~  *a  +  vi 

with  and  ~~  following  similarly. 
dyi  dz\ 


(A.41) 


Bearing  in  Remote  Deck  Coordinates 


For  this  final  instance,  where 


(A.26)  dictates  that 


x'2  =  Vii(a:i  +  A/?*i)  +  Vu(yi  +  Ai2ya) -f  ^3(^1  +  A/2,i) 
1/2  —  +  A/2,i)  +  Vaa(l/i  +  A/2„a)  +  Va3(zi  +  A/2,i) 

Thus,  the  derivatives  fake  the  form 


dB?  y'fVn+x'l'Vn 

d*x  «)2  +  (tf')a 


with  and  following  similarly. 
oyi  oz\ 


(A.42) 

(A.43) 

(A.44) 

(A.45) 


Tables  of  Mp  and  H 

Tables  A-l  through  A-4  formulas  cover  the  four  coordinate  systems  that  have  been 
introduced.  It  is  assumed  that  when  a  remote  coordinate  system  is  used  that 

A/2i  =  |  A/2*i  A/2yi  A/2,i  ]  (A.46) 

and  coordinate  rotation  matrix  N  have  been  precomputed.  Also,  when  a  deck  coordinate 
system  is  used,  precomputation  of  the  appropriate  coordinate  rotation  matrices 
and  H 'W  is  assumed. 

The  user  should  also  yote  that  these  formulations  involve  coordinate  transformations  of 
[  x  •••  y  •••  z  J  that  do  not  use  platform  velocity  as  input.  This  means  that  it 
is  implicitly  assumed  that  the  target  velocity  is  much  larger  than  that  of  the  ship.  While 
usually  true,  this  assumption  may  noticeably  degrade  the  otherwise  high  accuracy  of  doppler 
radar  measurements. 
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Table  A-l.  Formulas  for  Local  Stabilized  Coordinates 
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Table  A-2.  Formulas  for  Local  Deck  Coordinates 


Bearin 


s  -i  x'."  f 

M,  =  B["  tan  ^  H  = 
r  1  V?  L  dxx 

x'{'  ~  x\Un  4-  J/i U\2  4  *\U\3  yjw  =  4-  yit/jj  4-  sit^as 


9J5J"  A  as;"  afl',"  A 
0  0  ~d£~  0 


gflj"  ti'Un-xrUn  dB[" 

0*i  ■  (*i')a  4-  (yi")a  to  “  (*?')* +  (vi")J 


*iw  =  *i^si  4-  yidsa  4-  ril/aa 
tf'Uis-*?Un 

dxx  ~  (*H24-«)a 


Slant  range 


m, = «r = „ = v^f+iF+i?  h=[^„m:„m:0] 

g<,  _  4-  Vx'Un  4  £»£  _  x'l'Un  4-  yT^aa  4-  rf^as  &»'{'  _  4-  ygga;  4-  *("£/» 

5*1  «x  5yi  *i  dxx  «i 


ox  x  «x 


C)  Elevation 


My  =  E["  =  tan-1 

«  =  4-  (yi")3  4-  {xT? 

7=  vW  +  MT 


KT  +  W") 


H 

= 

dE["  _ 

\ 

dx  1 

a 

dE["  _ 

1 

ayi 

a 

1 

11 

a 

D)  Doppler 


M  -^-IX^+OTUM1  J5T- f  «*21  2fL  $£  M:  «C  £fT 

Mp  ~  "  dt  V{Xl  >  +  +  Ul  }  H  -  l  tox  dxx  dyx  dyx  dxx  dix 


x'x  —  *i U\\  +  iliUii  4-  xxUxs  if"  —  *1^21  4  (ixUix  4  iii/23 

-  _  ,,1112.111  .  -JH2JII  1  _ III  1111 

1  =  xi  xx  4  1/1  Vi  4  *1  *1 


da'"  _  xx  rjx  1 
dx  1  «i  "if” 

£1 

5ii  «i 


da'x  _  Vi  nyx 
dyx  ax  «f 

ds'l'  _  yx 

<9yi  *1 


<9ii  5yi  flyi  #*i  flii 
ij"  =  iit/31  4-  i)iUz2  4  iil/33 


da'"  z±  _  rjz^ 
dx  1  «i  V  sf 


da'"  _  xx 
dix  ~  ex 
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Table  A-3.  Formulas  for  Remote  Stabilized  Coordinates 
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Table  A-4.  Formulas  for  Remote  Deck  Coordinates 
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